/**
 * \copyright
 * Copyright (c) 2020, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 *
 */

/**************************************************************************
   MSHLib - Object:
   Task:
   Programing:
   08/2005 WW/OK Encapsulation from rf_ele_msh
   last modified
**************************************************************************/

#include "mathlib.h"
#include <cmath>
#include <cstdlib>  //WW
#include <float.h>  //WW
// MSHLib
// WW#include "MSHEnums.h" // KR 2010/11/15
#include "msh_elem.h"

namespace MeshLib
{
/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
CElem::CElem(size_t Index)
    : CCore(Index),
      normal_vector(NULL),
      /*geo_type(t), */ owner(NULL),
      ele_dim(1),
      nnodes(0),
      nnodesHQ(0),
      nodes(nnodes),
      nodes_index(nnodes),
      patch_index(0),
      transform_tensor(NULL)
{
    grid_adaptation = -1;
    volume = 0.0;
    face_index = -1;
    no_faces_on_surface = 0;
    gravity_center[0] = gravity_center[1] = gravity_center[2] = 0.0;
    normal_vector = NULL;
    area = 1.0;      // WW
    excavated = -1;  // WX

#if defined(USE_PETSC)  // || defined(using other parallel scheme). WW
    g_index = NULL;
#endif
}
/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
CElem::CElem() : CCore(0), normal_vector(NULL)
{
    selected = 0;
    // matgroup_view = 0;
    grid_adaptation = -1;
    nnodes = 0;
    nnodesHQ = 0;
    ele_dim = 1;  // Dimension of element
    patch_index = 0;
    //
    volume = 0.0;
    face_index = -1;
    no_faces_on_surface = 0;
    owner = NULL;
    nodes.resize(8);  // Nodes of face
    transform_tensor = NULL;
    gravity_center[0] = gravity_center[1] = gravity_center[2] = 0.0;
    area = 1.0;  // WW area = 1.0
    normal_vector = NULL;
    excavated = -1;  // WX

#if defined(USE_PETSC)  // || defined(using other parallel scheme). WW
    g_index = NULL;
#endif
}

/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
CElem::CElem(size_t Index, CElem* onwer, int Face)
    : CCore(Index), normal_vector(NULL), owner(onwer)
{
    int i, j, k, n, ne;
    int faceIndex_loc[10];
    int edgeIndex_loc[10] = {};
    no_faces_on_surface = 0;
    n = owner->GetElementFaceNodes(Face, faceIndex_loc);
    face_index = Face;
    gravity_center[0] = gravity_center[1] = gravity_center[2] = 0.0;
    transform_tensor = NULL;
    normal_vector = NULL;
    excavated = -1;  // WX
    //
    switch (owner->geo_type)
    {
        // case MshElemType::LINE:  // 1-D bar element //KR need not be
        // processed
        case MshElemType::QUAD:  // 2-D quadrilateral element
            this->setElementProperties(MshElemType::LINE);
            break;
        case MshElemType::HEXAHEDRON:  // 3-D hexahedral element
            this->setElementProperties(MshElemType::QUAD);
            break;
        case MshElemType::TRIANGLE:  // 2-D triagular element
            this->setElementProperties(MshElemType::LINE);
            break;
        case MshElemType::TETRAHEDRON:  // 3-D tetrahedral element
            this->setElementProperties(MshElemType::TRIANGLE);
            break;
        case MshElemType::PRISM:  // 3-D prismatic element
            if (Face < 2)         // top or bottom face of the prism
                this->setElementProperties(MshElemType::TRIANGLE);
            else  // side of the prism
                this->setElementProperties(MshElemType::QUAD);
            break;
        case MshElemType::PYRAMID:  // 3-D pyramid element
            if (Face < 1)           // bottom face
                this->setElementProperties(MshElemType::QUAD);
            else  // side faces
                this->setElementProperties(MshElemType::TRIANGLE);
            break;
        default:
            std::cerr << "CElem::CElem MshElemType not handled"
                      << "\n";
    }

    patch_index = owner->patch_index;
    quadratic = owner->quadratic;
    nodes_index.resize(n);
    nodes.resize(n);

    boundary_type = 'B';
    for (i = 0; i < n; i++)
    {
        nodes_index[i] = owner->nodes_index[faceIndex_loc[i]];
        nodes[i] = owner->nodes[faceIndex_loc[i]];
        // 18.02.2009. cf. changes in mapping & generation. WW
        if ((nodes[i]->GetBoundaryType() != '0') &&
            (nodes[i]->GetBoundaryType() != '1'))
            nodes[i]->SetBoundaryType('B');
    }
    // Face edges
    ne = owner->GetEdgesNumber();
    edges.resize(nnodes);
    edges_orientation.resize(nnodes);
    edges_orientation = 1;
    for (i = 0; i < nnodes; i++)
    {
        k = (i + 1) % nnodes;
        for (j = 0; j < ne; j++)
        {
            owner->GetLocalIndicesOfEdgeNodes(j, edgeIndex_loc);
            if ((faceIndex_loc[i] == edgeIndex_loc[0] &&
                 faceIndex_loc[k] == edgeIndex_loc[1]) ||
                (faceIndex_loc[i] == edgeIndex_loc[1] &&
                 faceIndex_loc[k] == edgeIndex_loc[0]))
            {
                edges[i] = owner->edges[j];
                if (faceIndex_loc[i] == edgeIndex_loc[1] &&
                    faceIndex_loc[k] == edgeIndex_loc[0])
                    edges_orientation[i] = -1;
                edges[i]->boundary_type = 'B';
                break;
            }
        }
    }

#if defined(USE_PETSC)  // || defined(using other parallel scheme). WW
    g_index = NULL;
#endif
}

/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW/OK Implementation
**************************************************************************/
CElem::CElem(size_t Index, CElem* m_ele_parent)
    : CCore(Index), normal_vector(NULL)
{
    //  static int faceIndex_loc[10];
    //  static int edgeIndex_loc[10];
    gravity_center[0] = gravity_center[1] = gravity_center[2] = 0.0;
    transform_tensor = NULL;
    this->setElementProperties(m_ele_parent->geo_type);

    patch_index = m_ele_parent->patch_index;
    quadratic = m_ele_parent->quadratic;

    boundary_type = 'I';

    // Initialize topological properties
    neighbors.resize(nfaces);
    for (size_t i = 0; i < nfaces; i++)
        neighbors[i] = NULL;
    edges.resize(nedges);
    edges_orientation.resize(nedges);
    for (size_t i = 0; i < nedges; i++)
    {
        edges[i] = NULL;
        edges_orientation[i] = 1;
    }
    area = 1.0;  // WW

    excavated = -1;  // 12.08.2011. WW

#if defined(USE_PETSC)  // || defined(using other parallel scheme). WW
    g_index = NULL;
#endif
}

CElem::CElem(CElem const& elem)
    : CCore(elem.GetIndex()),
      mat_vector(elem.mat_vector),
      // matgroup_view (elem.matgroup_view),
      selected(elem.selected),
      normal_vector(new double[3]),
      representative_length(elem.representative_length),
      courant(elem.courant),
      neumann(elem.neumann),
      geo_type(elem.geo_type),
      owner(elem.owner),
      ele_dim(elem.ele_dim),
      nnodes(elem.nnodes),
      nnodesHQ(elem.nnodesHQ),
      nodes(elem.nodes),
      nedges(elem.nedges),
      edges(elem.edges),
      edges_orientation(elem.edges_orientation),
      nfaces(elem.nfaces),
      no_faces_on_surface(elem.no_faces_on_surface),
      face_index(elem.face_index),
      volume(elem.volume),
      grid_adaptation(elem.grid_adaptation),
      patch_index(elem.patch_index),
      area(elem.area),
      angle(elem.angle)
{
    for (size_t k(0); k < 3; k++)
    {
        normal_vector[k] = elem.normal_vector[k];
        gravity_center[k] = elem.gravity_center[k];
    }

    // copy nodes
    nodes.resize((int)((elem.nodes).Size()));
    for (size_t k(0); k < (elem.nodes).Size(); k++)
        nodes[k] =
            new CNode((elem.nodes[k])->GetIndex(), (elem.nodes[k])->getData());

        // copy edges
        //	edges = vec<CEdge*> ((elem.edges).Size());
        //	for (size_t k(0); k<(elem.edges).Size() ; k++) {
        //		edges[k] = new CNode ((elem.nodes[k])->GetIndex(),
        //(elem.nodes[k])->X(), (elem.nodes[k])->Y(), (elem.nodes[k])->Z());
        //	}

        //*angle = *(elem.angle);
        // copy transform tensor
        //	Matrix * transform_tensor;

        // copy neighbors
        //	vec<CElem*> neighbors;

#if defined(USE_PETSC)  // || defined(using other parallel scheme). WW
    g_index = elem.g_index;
#endif
}

CElem::CElem(MshElemType::type t, size_t node0, size_t node1, size_t node2,
             int mat)
    : CCore(0),
      normal_vector(NULL),
      geo_type(t),
      owner(NULL),
      ele_dim(2),
      nnodes(3),
      nnodesHQ(6),
      nodes(nnodes),
      nodes_index(nnodes),
      nedges(3),
      edges(nedges),
      edges_orientation(nedges),
      nfaces(3),
      patch_index(mat),
      transform_tensor(NULL),
      neighbors(nfaces)
{
    nodes_index[0] = node0;
    nodes_index[1] = node1;
    nodes_index[2] = node2;

    // Initialize topological properties
    for (size_t i = 0; i < nfaces; i++)
        neighbors[i] = NULL;
    for (size_t i = 0; i < nedges; i++)
    {
        edges[i] = NULL;
        edges_orientation[i] = 1;
    }

#if defined(USE_PETSC)  // || defined(using other parallel scheme). WW
    g_index = NULL;
#endif
}

CElem::CElem(MshElemType::type t, size_t node0, size_t node1, size_t node2,
             size_t node3, int mat)
    : CCore(0),
      normal_vector(NULL),
      geo_type(t),
      owner(NULL),
      ele_dim(2),
      nnodes(4),
      nnodesHQ(9),
      nodes(nnodes),
      nodes_index(nnodes),
      nedges(4),
      edges(nedges),
      edges_orientation(nedges),
      nfaces(4),
      patch_index(mat),
      transform_tensor(NULL),
      neighbors(nfaces)
{
    nodes_index[0] = node0;
    nodes_index[1] = node1;
    nodes_index[2] = node2;
    nodes_index[3] = node3;

    // Initialize topological properties
    for (size_t i = 0; i < nfaces; i++)
        neighbors[i] = NULL;
    for (size_t i = 0; i < nedges; i++)
    {
        edges[i] = NULL;
        edges_orientation[i] = 1;
    }
#if defined(USE_PETSC)  // || defined(using other parallel scheme). WW
    g_index = NULL;
#endif
}

/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
CElem::~CElem()
{
    // HACK LB those resize() do not free memory allocated with new
    // HACK LB Possible memory leak otherwise resize() not necessary
    nodes_index.resize(0);
    nodes.resize(0);
    edges.resize(0);
    neighbors.resize(0);
    mat_vector.resize(0);
    edges_orientation.resize(0);
    owner = NULL;
    if (transform_tensor)
        delete transform_tensor;
    transform_tensor = NULL;
    transform_tensor = NULL;
    if (normal_vector)
        delete[] normal_vector;
    normal_vector = NULL;

#if defined(USE_PETSC)  // || defined(using other parallel scheme). WW
    if (g_index)
        delete[] g_index;
    g_index = NULL;
#endif
}
/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
void CElem::FillTransformMatrix(const bool recompute_matrix)
{
    if (transform_tensor && recompute_matrix == false)
        return;

    double xx[3];
    double yy[3];
    double zz[3];
    transform_tensor = new Math_Group::Matrix(3, 3);
    if (geo_type == MshElemType::LINE)
    {
        // x"_vec
        double const* const pnt0(nodes[0]->getData());
        double const* const pnt1(nodes[1]->getData());
        xx[0] = pnt1[0] - pnt0[0];
        xx[1] = pnt1[1] - pnt0[1];
        xx[2] = pnt1[2] - pnt0[2];
        NormalizeVector(xx, 3);
        // an arbitrary vector
        for (size_t i = 0; i < 3; i++)
            yy[i] = 0.0;
        // WW. 06.11.2007
        if (fabs(xx[0]) > 0.0 && fabs(xx[1]) + fabs(xx[2]) < DBL_MIN)
            yy[2] = 1.0;
        else if (fabs(xx[1]) > 0.0 && fabs(xx[0]) + fabs(xx[2]) < DBL_MIN)
            yy[0] = 1.0;
        else if (fabs(xx[2]) > 0.0 && fabs(xx[0]) + fabs(xx[1]) < DBL_MIN)
            yy[1] = 1.0;
        else
        {
            for (size_t i = 0; i < 3; i++)
                if (fabs(xx[i]) > 0.0)
                {
                    yy[i] = -xx[i];
                    break;
                }
        }
        // z"_vec
        CrossProduction(xx, yy, zz);
        NormalizeVector(zz, 3);
        // y"_vec
        CrossProduction(zz, xx, yy);
        NormalizeVector(yy, 3);
    }
    else if (geo_type == MshElemType::QUAD || geo_type == MshElemType::QUAD8 ||
             geo_type == MshElemType::TRIANGLE)
    {
        // x"_vec
        //			xx[0] = nodes[1]->X() - nodes[0]->X();
        //			xx[1] = nodes[1]->Y() - nodes[0]->Y();
        //			xx[2] = nodes[1]->Z() - nodes[0]->Z();
        double const* const pnt0(nodes[0]->getData());
        double const* const pnt1(nodes[1]->getData());
        xx[0] = pnt1[0] - pnt0[0];
        xx[1] = pnt1[1] - pnt0[1];
        xx[2] = pnt1[2] - pnt0[2];
        NormalizeVector(xx, 3);
        // a vector on the plane
        //			yy[0] = nodes[2]->X() - nodes[1]->X();
        //			yy[1] = nodes[2]->Y() - nodes[1]->Y();
        //			yy[2] = nodes[2]->Z() - nodes[1]->Z();
        double const* const pnt2(nodes[2]->getData());
        yy[0] = pnt2[0] - pnt1[0];
        yy[1] = pnt2[1] - pnt1[1];
        yy[2] = pnt2[2] - pnt1[2];
        // z"_vec. off plane
        CrossProduction(xx, yy, zz);
        NormalizeVector(zz, 3);
        // y"_vec
        CrossProduction(zz, xx, yy);
        NormalizeVector(yy, 3);
    }
    for (size_t i = 0; i < 3; i++)
    {
        (*transform_tensor)(i, 0) = xx[i];
        (*transform_tensor)(i, 1) = yy[i];
        (*transform_tensor)(i, 2) = zz[i];
    }
}
/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   02/2006 PCH Implementation
**************************************************************************/
double CElem::getTransformTensor(int idx)
{
    // WW
    int i = idx % 3;
    int j = idx / 3;
    return (*transform_tensor)(i, j);
    // return MatT[idx];
}

/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
void CElem::SetFace()
{
    nodes.resize(8);
    nodes_index.resize(8);
}

/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
void CElem::SetFace(CElem* onwer, const int Face)
{
    static int nodeIndex_loc[8];
    no_faces_on_surface = 0;
    owner = onwer;
    size_t n = owner->GetElementFaceNodes(Face, nodeIndex_loc);
    quadratic = owner->quadratic;
    face_index = Face;
    patch_index = owner->patch_index;
    switch (owner->geo_type)
    {
        // case MshElemType::LINE:  // 1-D bar element
        case MshElemType::QUAD:  // 2-D quadrilateral element
            this->setElementProperties(MshElemType::LINE);  // JOD 2014-11-10
            break;
        case MshElemType::HEXAHEDRON:  // 3-D hexahedral element
            this->setElementProperties(MshElemType::QUAD8);
            break;
        // case MshElemType::TRIANGLE:  // 2-D triagular element
        case MshElemType::TETRAHEDRON:  // 3-D tetrahedral element
            this->setElementProperties(MshElemType::TRIANGLE);
            break;
        case MshElemType::PRISM:
            if (Face < 2)
                this->setElementProperties(MshElemType::TRIANGLE);
            else
                this->setElementProperties(MshElemType::QUAD8);
            break;  // 3-D prismatic element
        case MshElemType::PYRAMID:
            if (Face < 1)
                this->setElementProperties(MshElemType::QUAD8);
            else
                this->setElementProperties(MshElemType::TRIANGLE);
            break;  // 3-D pyramid element
        default:
            std::cerr << "CElem::SetFace MshElemType not handled"
                      << "\n";
            break;
    }

    if (nodes.Size() < n)
        nodes.resize(n);
    for (size_t i = 0; i < n; i++)
    {
        nodes[i] = owner->nodes[nodeIndex_loc[i]];
        nodes_index[i] = nodes[i]->GetIndex();
    }
}
/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
   11/2010 KR moved to MSHEnums.h
**************************************************************************/
std::string CElem::GetName() const
{
    return MshElemType2String(geo_type);
}

/**************************************************************************
   MSHLib-Method:
   Programing:
   06/2005 WW Implementation
   08/2005 WW/OK Extension for GMS/SOL files
   10/2008 OK FEFLOW
**************************************************************************/
void CElem::Read(std::istream& is, int fileType)
{
    // fileType=0: msh
    // fileType=1: rfi
    // fileType=2: gmsh
    // fileType=3: GMS
    // fileType=4: SOL
    // fileType=5: FLAC 3D. WW
    // fileType=5x: FLAC3D //MR
    // fileType=6: FEFLOW //OK
    // fileType=7: GMSH 2008 Version//TK
    int idummy(-1), et(-1);
    std::string name("");
    int gmsh_patch_index;  // TKOK

    //   is.ignore(numeric_limits<int>::max(), '\n');
    //----------------------------------------------------------------------
    // 1 Reading element type data
    switch (fileType)
    {
        //....................................................................
        case 0:  // msh
        {
            std::string buffer("");
            is >> index >> patch_index;
            is >> buffer;

            if (buffer.find("-1") != std::string::npos)
            {
                grid_adaptation = strtol(buffer.data(), NULL, 0);
                is >> name;
            }
            else
                name = buffer;

            geo_type = String2MshElemType(name);
            break;
        }
        //....................................................................
        case 1:  // rfi
            is >> index >> patch_index >> name;
            geo_type = String2MshElemType(name);
            break;
        //....................................................................
        case 2:  // gmsh
            is >> index >> et >> gmsh_patch_index >> idummy >> nnodes;
            patch_index = gmsh_patch_index - 1;  // OK
            switch (et)
            {
                case 1:
                    geo_type = MshElemType::LINE;
                    break;
                case 2:
                    geo_type = MshElemType::TRIANGLE;
                    break;
                case 3:
                    geo_type = MshElemType::QUAD;
                    break;
                case 4:
                    geo_type = MshElemType::TETRAHEDRON;
                    break;
                case 5:
                    geo_type = MshElemType::HEXAHEDRON;
                    break;
                case 6:
                    geo_type = MshElemType::PRISM;
                    break;
                case 7:
                    geo_type = MshElemType::PYRAMID;
                    break;
            }
            index--;
            break;
        case 7:  // GMSH 2008
            size_t nb_tags;

            is >> index >> et >> nb_tags >> idummy >> gmsh_patch_index;
            patch_index = gmsh_patch_index;
            for (size_t j = 2; j < nb_tags; j++)
                is >> idummy;
            switch (et)
            {
                case 1:
                    geo_type = MshElemType::LINE;
                    patch_index =
                        0;  // KR: can line elements have material ids?
                    nnodes = 2;
                    break;
                case 2:
                    geo_type = MshElemType::TRIANGLE;
                    nnodes = 3;
                    break;
                case 3:
                    geo_type = MshElemType::QUAD;
                    nnodes = 4;
                    break;
                case 4:
                    geo_type = MshElemType::TETRAHEDRON;
                    nnodes = 4;
                    break;
                case 5:
                    geo_type = MshElemType::HEXAHEDRON;
                    nnodes = 8;
                    break;
                case 6:
                    geo_type = MshElemType::PRISM;
                    nnodes = 6;
                    break;
                case 7:
                    geo_type = MshElemType::PYRAMID;
                    break;
                case 15:
                    geo_type = MshElemType::INVALID;
                    nnodes = 1;
                    break;
                default:
                    geo_type = MshElemType::INVALID;
                    break;
            }
            index--;
            break;
        //....................................................................
        case 3:  // GMS
            geo_type = MshElemType::TRIANGLE;
            break;
        //....................................................................
        case 4:  // gmsh
            geo_type = MshElemType::TRIANGLE;
            break;
        //....................................................................
        case 5:  // FLAC 3D. 14.01.2008 WW
            geo_type = MshElemType::HEXAHEDRON;
            break;
        case 56:                            // FLAC3D - pri (Wedge)        //MR
            geo_type = MshElemType::PRISM;  // MR
            fileType = 5;                   // MR
            break;                          // MR
    }

    // TF
    if (geo_type == MshElemType::INVALID)
    {
        // read rest of line
        std::string tmp;
        getline(is, tmp);
        return;
    }

    //----------------------------------------------------------------------
    // 2 Element configuration
    this->setElementProperties(geo_type);

    //----------------------------------------------------------------------
    // 3 Reading element node data
    switch (fileType)
    {
        case 0:  // msh
            for (int i = 0; i < nnodes; i++)
                is >> nodes_index[i];
            break;
        case 1:  // rfi
            for (int i = 0; i < nnodes; i++)
                is >> nodes_index[i];
            break;
        case 2:  // gmsh
            for (int i = 0; i < nnodes; i++)
            {
                is >> nodes_index[i];
                nodes_index[i] -= 1;
            }
            break;
        case 7:  // GMSH 2008
            if (et != 15)
                for (int i = 0; i < nnodes; i++)
                {
                    is >> nodes_index[i];
                    nodes_index[i] -= 1;
                }
            else
            {
                // eat rest of line
                std::string dummy;
                is >> dummy;
            }
            break;
        case 3:  // GMS
            for (int i = 0; i < nnodes; i++)
            {
                is >> nodes_index[i];
                nodes_index[i] -= 1;
            }
            break;
        case 4:  // SOL
            for (int i = 0; i < nnodes; i++)
            {
                is >> nodes_index[i];
                nodes_index[i] -= 1;
            }
            is >> patch_index;
            break;
        case 5:  // FLAC 3D. 14.01.2008. WW
            for (int i = 0; i < nnodes; i++)
            {
                is >> nodes_index[i];
                nodes_index[i] -= 1;
            }
            break;
        case 6:  // FEFLOLW
            for (int i = 0; i < nnodes; i++)
            {
                is >> nodes_index[i];
                nodes_index[i] -= 1;
            }
            break;
        case 8:  // GMS_3DM
            is >> idummy;
            for (int i = 0; i < nnodes; i++)
            {
                is >> nodes_index[i];
                nodes_index[i] -= 1;
            }
            is >> patch_index;
            patch_index -= 1;
            break;
    }
    is >> std::ws;
    //----------------------------------------------------------------------
    // Initialize topological properties
    neighbors.resize(nfaces);
    for (size_t i = 0; i < nfaces; i++)
        neighbors[i] = NULL;
    edges.resize(nedges);
    edges_orientation.resize(nedges);
    for (size_t i = 0; i < nedges; i++)
    {
        edges[i] = NULL;
        edges_orientation[i] = 1;
    }
}

void CElem::InitializeMembers()
{
    // Initialize topological properties
    neighbors.resize(nfaces);
    for (size_t i = 0; i < nfaces; i++)
        neighbors[i] = NULL;
    edges.resize(nedges);
    edges_orientation.resize(nedges);
    for (size_t i = 0; i < nedges; i++)
    {
        edges[i] = NULL;
        edges_orientation[i] = 1;
    }
}

/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
void CElem::WriteIndex(std::ostream& os) const
{
    // Comment for GUI WW if(quadratic) nn = nnodesHQ;
    os << index << " " << patch_index << " " << GetName() << " ";
    for (int i = 0; i < nnodes - 1; i++)
        os << nodes_index[i] << " ";
    os << nodes_index[nnodes - 1] << "\n";
}

/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
void CElem::WriteIndex_TEC(std::ostream& os) const
{
    std::string deli = "  ";
    if (geo_type == MshElemType::LINE)
        os << nodes_index[0] + 1 << deli << nodes_index[1] + 1 << deli
           << nodes_index[1] + 1 << deli << nodes_index[0] + 1;

    else if (geo_type == MshElemType::TRIANGLE)
        os << nodes_index[0] + 1 << deli << nodes_index[1] + 1 << deli
           << nodes_index[2] + 1 << deli << nodes_index[0] + 1;
    else if (geo_type == MshElemType::PRISM)
        os << nodes_index[0] + 1 << deli << nodes_index[0] + 1 << deli
           << nodes_index[1] + 1 << deli << nodes_index[2] + 1 << deli
           << nodes_index[3] + 1 << deli << nodes_index[3] + 1 << deli
           << nodes_index[4] + 1 << deli << nodes_index[5] + 1 << deli;
    else if (geo_type == MshElemType::PYRAMID)
        os << nodes_index[0] + 1 << deli << nodes_index[1] + 1 << deli
           << nodes_index[2] + 1 << deli << nodes_index[3] + 1 << deli
           << nodes_index[4] + 1 << deli << nodes_index[4] + 1 << deli
           << nodes_index[4] + 1 << deli << nodes_index[4] + 1 << deli;
    else

        for (int i = 0; i < nnodes; i++)
            os << nodes_index[i] + 1 << deli;
    os << '\n';  // GK44 for io performance do not flush buffer with endl...
}
/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
void CElem::WriteAll(std::ostream& os) const
{
    std::string deli = "  ";
    os << index << deli << patch_index << deli << GetName() << deli;
    os << "Index X Y Z: "
       << "\n";
    for (size_t i = 0; i < nodes.Size(); i++)
    {
        double const* const pnt_i(nodes[i]->getData());
        os << nodes_index[i] << deli << pnt_i[0] << deli << pnt_i[1] << deli
           << pnt_i[2] << "\n";
    }
}

/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
void CElem::WriteNeighbors(std::ostream& os) const
{
    os << "Neighbors of " << index << "\n";
    for (size_t i = 0; i < nfaces; i++)
        neighbors[i]->WriteAll(os);
    os << "End neighbors of " << index << "\n"
       << "\n";
}

/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
void CElem::MarkingAll(bool makop)
{
    this->mark = makop;
    int SizeV = nnodes;
    if (quadratic)
        SizeV = nnodesHQ;
    for (int i = 0; i < SizeV; i++)
        nodes[i]->SetMark(makop);

    size_t nedg = edges.Size();
    for (size_t i = 0; i < nedg; i++)
        edges[i]->SetMark(makop);
}
/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
void CElem::SetNodes(Math_Group::vec<CNode*>& ele_nodes, bool ReSize)
{
    int SizeV = nnodes;
    if (quadratic)
        SizeV = nnodesHQ;
    if (ReSize)
    {
        nodes.resize(SizeV);
#if !defined(USE_PETSC)  // && !defined (other parallel linear solver lib).
                         // //WW. 05.2013
        nodes_index.resize(SizeV);
#endif
    }
    for (int i = 0; i < SizeV; i++)
    {
        nodes[i] = ele_nodes[i];
#if !defined(USE_PETSC)  // && !defined (other parallel linear solver lib).
                         // //WW. 05.2013
        nodes_index[i] = nodes[i]->GetIndex();
#endif
    }
}

void CElem::setNodes(std::vector<CNode*> const& ele_nodes)
{
    if (nodes.Size() != ele_nodes.size())
    {
        nodes.resize(ele_nodes.size());
        nodes_index.resize(ele_nodes.size());
    }

    for (size_t k(0); k < ele_nodes.size(); k++)
    {
        nodes[k] = ele_nodes[k];
        nodes_index[k] = ele_nodes[k]->GetIndex();
    }
}

/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   06/2005 WW Implementation
**************************************************************************/
void CElem::GetLocalIndicesOfEdgeNodes(const int Edge, int* EdgeNodes)
{
    switch (geo_type)
    {
        case MshElemType::LINE:
            EdgeNodes[0] = 0;
            EdgeNodes[1] = 1;
            break;               // 1-D bar element
        case MshElemType::QUAD:  // 2-D quadrilateral element
            EdgeNodes[0] = Edge;
            EdgeNodes[1] = (Edge + 1) % 4;
            break;
        case MshElemType::HEXAHEDRON:  // 3-D hexahedral element
            if (Edge < 8)
            {
                EdgeNodes[0] = Edge;
                EdgeNodes[1] = (Edge + 1) % 4 + 4 * (int)(Edge / 4);
            }
            else
            {
                EdgeNodes[0] = Edge % 4;
                EdgeNodes[1] = Edge % 4 + 4;
            }
            break;
        case MshElemType::TRIANGLE:  // 2-D triagular element
            EdgeNodes[0] = Edge;
            EdgeNodes[1] = (Edge + 1) % 3;
            break;
        case MshElemType::TETRAHEDRON:  // 3-D tetrahedra
            if (Edge < 3)
            {
                EdgeNodes[0] = Edge;
                EdgeNodes[1] = (Edge + 1) % 3;
            }
            else
            {
                EdgeNodes[0] = 3;
                EdgeNodes[1] = (Edge + 1) % 3;
            }

            break;
        case MshElemType::PRISM:  // 3-D prismatic element
            if (Edge < 6)
            {
                EdgeNodes[0] = Edge;
                EdgeNodes[1] = (Edge + 1) % 3 + 3 * (int)(Edge / 3);
            }
            else
            {
                EdgeNodes[0] = Edge % 3;
                EdgeNodes[1] = Edge % 3 + 3;
            }
            break;
        case MshElemType::PYRAMID:  // 3-D pyramid element
            if (Edge < 4)
            {
                EdgeNodes[0] = Edge;
                EdgeNodes[1] = (Edge + 1) % 4;
            }
            else
            {
                EdgeNodes[0] = Edge % 4;
                EdgeNodes[1] = 4;
            }
            break;
        default:
            std::cerr << "CElem::GetLocalIndicesOfEdgeNodes() - MshElemType "
                         "not handled"
                      << "\n";
    }
}
/**************************************************************************
   GetElementFaceNodes
   Task: Get local indeces of an element face nodes
   Return: number of nodes of a face
   Programing:
   06/2004 WW
**************************************************************************/
int CElem::GetElementFaces1D(int* FaceNode)
{
    FaceNode[0] = 0;
    FaceNode[1] = 1;
    return 2;
}
/**************************************************************************
   GetElementFaceNodesTri
   Task: Get local indeces of a traingle element face nodes
   Augs.:
        const int Face :  Local index of element face
        const int order:  1 Linear. 2, quadratic
        int *FaceNode  :  Local index of face nodes
   Return: number of nodes of a face
   Programing:
   09/2004 WW
**************************************************************************/
int CElem::GetElementFacesTri(int Face, int* FaceNode)
{
    if (!quadratic)
    {
        FaceNode[0] = Face;
        FaceNode[1] = (Face + 1) % 3;
        return 2;
    }
    else
    {
        FaceNode[0] = Face;
        FaceNode[1] = (Face + 1) % 3;
        FaceNode[2] = Face + 3;
        return 3;
    }
}
/**************************************************************************
   GetElementFaceNodesQuad
   Task: Get local indeces of a quadralateral element face nodes
   Augs.:
        const int Face :  Local index of element face
        const int order:  1 Linear. 2, quadratic
        int *FaceNode  :  Local index of face nodes
   Return: number of nodes of a face
   Programing:
   09/2004 WW
**************************************************************************/
int CElem::GetElementFacesQuad(int Face, int* FaceNode)
{
    if (!quadratic)
    {
        FaceNode[0] = Face;
        FaceNode[1] = (Face + 1) % 4;
        return 2;
    }
    else
    {
        FaceNode[0] = Face;
        FaceNode[1] = (Face + 1) % 4;
        FaceNode[2] = Face + 4;
        return 3;
    }
}
/**************************************************************************
   GetElementFaceNodesHex
   Task: Get local indeces of a hexahedra element face nodes
   Augs.:
        const int Face :  Local index of element face
        const int order:  1 Linear. 2, quadratic
        int *FaceNode  :  Local index of face nodes
   Return: number of nodes of a face
   Programing:
   09/2004 WW
**************************************************************************/
int CElem::GetElementFacesHex(int Face, int* FaceNode)
{
    int nn = 4, k = 0;
    if (quadratic)
        nn = 8;
    switch (Face)
    {
        case 0:
            for (k = 0; k < 4; k++)
                FaceNode[k] = k;
            if (quadratic)
                for (k = 0; k < 4; k++)
                    FaceNode[k + 4] = k + 8;
            break;
        case 1:
            for (k = 0; k < 4; k++)
                FaceNode[k] = k + 4;
            if (quadratic)
                for (k = 0; k < 4; k++)
                    FaceNode[k + 4] = k + 12;
            break;
        case 2:
            FaceNode[0] = 0;
            FaceNode[1] = 4;
            FaceNode[2] = 5;
            FaceNode[3] = 1;
            if (quadratic)
            {
                FaceNode[4] = 16;
                FaceNode[5] = 12;
                FaceNode[6] = 17;
                FaceNode[7] = 8;
            }
            break;
        case 3:
            FaceNode[0] = 1;
            FaceNode[1] = 5;
            FaceNode[2] = 6;
            FaceNode[3] = 2;
            if (quadratic)
            {
                FaceNode[4] = 17;
                FaceNode[5] = 13;
                FaceNode[6] = 18;
                FaceNode[7] = 9;
            }

            break;
        case 4:
            FaceNode[0] = 2;
            FaceNode[1] = 6;
            FaceNode[2] = 7;
            FaceNode[3] = 3;
            if (quadratic)
            {
                FaceNode[4] = 18;
                FaceNode[5] = 14;
                FaceNode[6] = 19;
                FaceNode[7] = 10;
            }
            break;
        case 5:
            FaceNode[0] = 0;
            FaceNode[1] = 3;
            FaceNode[2] = 7;
            FaceNode[3] = 4;
            if (quadratic)
            {
                FaceNode[4] = 11;
                FaceNode[5] = 19;
                FaceNode[6] = 15;
                FaceNode[7] = 16;
            }
            break;
    }
    return nn;
}
/**************************************************************************
   GetElementFaceNodesTet
   Task: Get local indeces of a Tedrahedra element face nodes
   Augs.:
        const int Face :  Local index of element face
        const int order:  1 Linear. 2, quadratic
        int *FaceNode  :  Local index of face nodes
   Return: number of nodes of a face
   Programing:
   09/2004 WW
**************************************************************************/
int CElem::GetElementFacesTet(int Face, int* FaceNode)
{
    int nn = 3;
    if (quadratic)
        nn = 6;
    switch (Face)
    {
        case 0:
            FaceNode[0] = 1;
            FaceNode[1] = 2;
            FaceNode[2] = 3;
            if (quadratic)
            {
                FaceNode[3] = 5;
                FaceNode[4] = 8;
                FaceNode[5] = 7;
            }
            break;
        case 1:
            FaceNode[0] = 3;
            FaceNode[1] = 2;
            FaceNode[2] = 0;
            if (quadratic)
            {
                FaceNode[3] = 8;
                FaceNode[4] = 6;
                FaceNode[5] = 9;
            }
            break;
        case 2:
            FaceNode[0] = 1;
            FaceNode[1] = 3;
            FaceNode[2] = 0;
            if (quadratic)
            {
                FaceNode[3] = 7;
                FaceNode[4] = 9;
                FaceNode[5] = 4;
            }
            break;
        case 3:
            FaceNode[0] = 0;
            FaceNode[1] = 2;
            FaceNode[2] = 1;
            if (quadratic)
            {
                FaceNode[3] = 6;
                FaceNode[4] = 5;
                FaceNode[5] = 4;
            }
            break;
    }
    return nn;
}

/**************************************************************************
   GetElementFaceNodesPri
   Task: Get local indeces of a prismal element face nodes
   Augs.:
        int Face :  Local index of element face
        int *FaceNode  :  Local index of face nodes
   Return: number of nodes of a face
   Programing:
   09/2004 WW
**************************************************************************/
int CElem::GetElementFacesPri(int Face, int* FaceNode)
{
    int nn = 3, k = 0;
    switch (Face)
    {
        case 0:
            nn = 3;
            for (k = 0; k < 3; k++)
                FaceNode[k] = k;
            if (quadratic)
            {
                for (k = 0; k < 3; k++)
                    FaceNode[k + 3] = k + 6;
                nn = 6;
            }
            break;
        case 1:
            for (k = 0; k < 3; k++)
                FaceNode[k] = k + 3;
            nn = 3;
            if (quadratic)
            {
                for (k = 0; k < 3; k++)
                    FaceNode[k + 3] = k + 9;
                nn = 6;
            }
            break;
        case 2:
            FaceNode[0] = 1;
            FaceNode[1] = 2;
            FaceNode[2] = 5;
            FaceNode[3] = 4;
            nn = 4;
            if (quadratic)
            {
                FaceNode[4] = 7;
                FaceNode[5] = 14;
                FaceNode[6] = 10;
                FaceNode[7] = 13;
                nn = 8;
            }
            break;
        case 3:
            FaceNode[0] = 5;
            FaceNode[1] = 2;
            FaceNode[2] = 0;
            FaceNode[3] = 3;
            nn = 4;
            if (quadratic)
            {
                FaceNode[4] = 14;
                FaceNode[5] = 8;
                FaceNode[6] = 12;
                FaceNode[7] = 11;
                nn = 8;
            }
            break;
        case 4:
            FaceNode[0] = 0;
            FaceNode[1] = 1;
            FaceNode[2] = 4;
            FaceNode[3] = 3;
            nn = 4;
            if (quadratic)
            {
                FaceNode[4] = 6;
                FaceNode[5] = 13;
                FaceNode[6] = 9;
                FaceNode[7] = 12;
                nn = 8;
            }
            break;
    }
    return nn;
}
/**************************************************************************
   GetElementFacesPyra
   Task: Get local indeces of a element face nodes
   Augs.:
        const int Face :  Local index of element face
        const int order:  1 Linear. 2, quadratic
        int *FaceNode  :  Local index of face nodes
   Return: number of nodes of a face
   Programing:
   01/2010 NW
**************************************************************************/
int CElem::GetElementFacesPyramid(const int Face, int* FaceNode)
{
    int nn = 3, k = 0;
    switch (Face)
    {
        case 0:
            nn = 4;
            for (k = 0; k < nn; k++)
                FaceNode[k] = k;
            if (quadratic)
            {
                for (k = 0; k < nn; k++)
                    FaceNode[k + 4] = k + 5;
                nn = 8;
            }
            break;
        case 1:
            FaceNode[0] = 0;
            FaceNode[1] = 1;
            FaceNode[2] = 4;
            nn = 3;
            if (quadratic)
            {
                FaceNode[3] = 5;
                FaceNode[4] = 10;
                FaceNode[5] = 9;
                nn = 6;
            }
            break;
        case 2:
            FaceNode[0] = 1;
            FaceNode[1] = 2;
            FaceNode[2] = 4;
            nn = 3;
            if (quadratic)
            {
                FaceNode[3] = 6;
                FaceNode[4] = 11;
                FaceNode[5] = 10;
                nn = 6;
            }
            break;
        case 3:
            FaceNode[0] = 2;
            FaceNode[1] = 3;
            FaceNode[2] = 4;
            nn = 3;
            if (quadratic)
            {
                FaceNode[3] = 7;
                FaceNode[4] = 12;
                FaceNode[5] = 11;
                nn = 6;
            }
            break;
        case 4:
            FaceNode[0] = 3;
            FaceNode[1] = 0;
            FaceNode[2] = 4;
            nn = 3;
            if (quadratic)
            {
                FaceNode[3] = 8;
                FaceNode[4] = 9;
                FaceNode[5] = 12;
                nn = 6;
            }
            break;
    }
    return nn;
}
/**************************************************************************
   GetElementFaces
   Task: set element faces (Geometry)
   Augs.:
        const int Face :  Local index of element face
        const int order:  1 Linear. 2, quadratic
        int *FaceNode  :  Local index of face nodes

   Programing:
   09/2004 WW
**************************************************************************/
int CElem::GetElementFaceNodes(int Face, int* FacesNode)
{
    switch (geo_type)
    {
        case MshElemType::LINE:  // 1-D bar element
            return GetElementFaces1D(FacesNode);
        case MshElemType::QUAD:  // 2-D quadrilateral element
            return GetElementFacesQuad(Face, FacesNode);
        case MshElemType::HEXAHEDRON:  // 3-D hexahedral element
            return GetElementFacesHex(Face, FacesNode);
        case MshElemType::TRIANGLE:  // 2-D triagular element
            return GetElementFacesTri(Face, FacesNode);
        case MshElemType::TETRAHEDRON:  // 3-D tetrahedral element
            return GetElementFacesTet(Face, FacesNode);
        case MshElemType::PRISM:
            return GetElementFacesPri(Face, FacesNode);
        // 3-D prismatic element
        case MshElemType::PYRAMID:  // 3-D pyramid element
            return GetElementFacesPyramid(Face, FacesNode);
        default:
            std::cerr << "CElem::GetElementFaceNodes MshElemType not handled"
                      << "\n";
    }
    return 0;
}

/**************************************************************************
   MSHLib-Method:
   Task: Compute volume for elements with straight edges and surfaces
   Programing:
   06/2005 WW Implementation
   03/2011 KR cleaned up code
**************************************************************************/
void CElem::ComputeVolume()
{
    volume = calcVolume();

    if (this->geo_type == MshElemType::LINE)  // Line
        representative_length = volume;
    else if (this->geo_type == MshElemType::TRIANGLE)
        representative_length = sqrt(volume) * 4.0;
    else if (this->geo_type == MshElemType::QUAD)
        representative_length = sqrt(volume);  // kg44 reactivated
    else if (this->geo_type == MshElemType::TETRAHEDRON)
        representative_length = sqrt(volume) * 6.0;
    else if (this->geo_type == MshElemType::HEXAHEDRON)
        representative_length = pow(volume, 1. / 3.);
    else if (this->geo_type == MshElemType::PRISM)
        representative_length = pow(volume, 1. / 3.);
    else if (this->geo_type == MshElemType::PYRAMID)
        representative_length = 0.0;  // NW set zero because I don't know what
                                      // is the representative_length.
    else
        std::cerr << "Error in CElem::ComputeVolume() - MshElemType not found"
                  << "\n";
}

double CElem::calcVolume() const
{
    double elemVolume = 0.0;

    if (this->geo_type == MshElemType::LINE)  // Line
    {
        double const* const pnt0(nodes[0]->getData());
        double const* const pnt(nodes[nnodes - 1]->getData());
        double xDiff = pnt[0] - pnt0[0];
        double yDiff = pnt[1] - pnt0[1];
        double zDiff = pnt[2] - pnt0[2];
        elemVolume = sqrt(xDiff * xDiff + yDiff * yDiff +
                          zDiff * zDiff);  // CMCD kg44 reactivated
    }
    else if (this->geo_type == MshElemType::TRIANGLE)
        elemVolume = ComputeDetTri(nodes[0]->getData(), nodes[1]->getData(),
                                   nodes[2]->getData());  // kg44 reactivated
    else if (this->geo_type == MshElemType::QUAD)
        elemVolume = ComputeDetTri(nodes[0]->getData(), nodes[1]->getData(),
                                   nodes[2]->getData()) +
                     ComputeDetTri(nodes[2]->getData(), nodes[3]->getData(),
                                   nodes[0]->getData());
    else if (this->geo_type == MshElemType::TETRAHEDRON)
        elemVolume = ComputeDetTex(nodes[0]->getData(), nodes[1]->getData(),
                                   nodes[2]->getData(),
                                   nodes[3]->getData());  // kg44 reactivated
    else if (this->geo_type == MshElemType::HEXAHEDRON)
    {
        elemVolume = ComputeDetTex(nodes[4]->getData(), nodes[7]->getData(),
                                   nodes[5]->getData(), nodes[0]->getData());
        elemVolume += ComputeDetTex(nodes[5]->getData(), nodes[3]->getData(),
                                    nodes[1]->getData(), nodes[0]->getData());
        elemVolume += ComputeDetTex(nodes[5]->getData(), nodes[7]->getData(),
                                    nodes[3]->getData(), nodes[0]->getData());
        elemVolume += ComputeDetTex(nodes[5]->getData(), nodes[7]->getData(),
                                    nodes[6]->getData(), nodes[2]->getData());
        elemVolume += ComputeDetTex(nodes[1]->getData(), nodes[3]->getData(),
                                    nodes[5]->getData(), nodes[2]->getData());
        elemVolume += ComputeDetTex(nodes[3]->getData(), nodes[7]->getData(),
                                    nodes[5]->getData(), nodes[2]->getData());
        // kg44 reactivated
    }
    else if (this->geo_type == MshElemType::PRISM)
    {
        elemVolume = ComputeDetTex(nodes[0]->getData(), nodes[1]->getData(),
                                   nodes[2]->getData(), nodes[3]->getData());
        elemVolume += ComputeDetTex(nodes[1]->getData(), nodes[4]->getData(),
                                    nodes[2]->getData(), nodes[3]->getData());
        elemVolume += ComputeDetTex(nodes[2]->getData(), nodes[4]->getData(),
                                    nodes[5]->getData(), nodes[3]->getData());
        // kg44 reactivated ---------Here the direction of flow needs to be
        // taken into account, we need rep length in x,y,z direction
    }
    else if (this->geo_type == MshElemType::PYRAMID)
    {
        elemVolume = ComputeDetTex(nodes[0]->getData(), nodes[1]->getData(),
                                   nodes[2]->getData(), nodes[4]->getData());
        elemVolume += ComputeDetTex(nodes[2]->getData(), nodes[3]->getData(),
                                    nodes[0]->getData(), nodes[4]->getData());
    }
    else
        std::cerr << "Error in CElem::ComputeVolume() - MshElemType not found"
                  << "\n";

    return elemVolume;
}

/**************************************************************************
   MSHLib-Method:
   Task:
   Programing:
   01/2006 YD Implementation
**************************************************************************/
void CElem::FaceNormal(int index0, int index1, double* face)
{
    if (GetElementType() == MshElemType::TRIANGLE ||
        GetElementType() == MshElemType::QUAD)
    {
        double xx[3];
        double yy[3];
        double zz[3];

        //----plane normal----------------------------
        // transform_tensor = new Matrix(3,3);
        // face"_vec
        double const* pnt0(nodes[index0]->getData());
        double const* pnt1(nodes[index1]->getData());
        face[0] = pnt1[0] - pnt0[0];
        face[1] = pnt1[1] - pnt0[1];
        face[2] = pnt1[2] - pnt0[2];

        // x"_vec
        //         xx[0] = nodes[1]->X()-nodes[0]->X();
        //         xx[1] = nodes[1]->Y()-nodes[0]->Y();
        //         xx[2] = nodes[1]->Z()-nodes[0]->Z();
        pnt0 = nodes[0]->getData();
        pnt1 = nodes[1]->getData();
        xx[0] = pnt1[0] - pnt0[0];
        xx[1] = pnt1[1] - pnt0[1];
        xx[2] = pnt1[2] - pnt0[2];
        // NormalizeVector(xx,3);
        // a vector on the plane
        //         yy[0] = nodes[2]->X()-nodes[1]->X();
        //         yy[1] = nodes[2]->Y()-nodes[1]->Y();
        //         yy[2] = nodes[2]->Z()-nodes[1]->Z();
        double const* pnt2(nodes[2]->getData());
        yy[0] = pnt2[0] - pnt1[0];
        yy[1] = pnt2[1] - pnt1[1];
        yy[2] = pnt2[2] - pnt1[2];
        // z"_vec

        CrossProduction(xx, yy, zz);
        NormalizeVector(zz, 3);
        // y"_vec
        CrossProduction(face, zz, yy);
        NormalizeVector(yy, 3);
        for (size_t i = 0; i < 3; i++)
            face[i] = yy[i];
    }
}

/**************************************************************************
   MSHLib-Method:
   06/2006 OK Implementation
**************************************************************************/
void CElem::SetNormalVector()
{
    if (!normal_vector)
        normal_vector = new double[3];                // WW
    if (this->GetElementType() == MshElemType::LINE)  // JOD 2014-11-10
    {
        double const* const p0(nodes[0]->getData());
        double const* const p1(nodes[1]->getData());
        double v1[3] = {p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]};
        if (fabs(v1[2]) > 1.e-20)
        {
            double buffer = v1[1];
            v1[1] = v1[2];
            v1[2] = buffer;
        }
        const double v2[3] = {0, 0, 1};  // fluxes are on xy plane
        CrossProduction(v1, v2, normal_vector);
        NormalizeVector(normal_vector, 3);
        if (normal_vector[0] < 0)
            normal_vector[0] = -normal_vector[0];
    }
    if (this->GetElementType() == MshElemType::TRIANGLE)
    {
        double const* const p0(nodes[0]->getData());
        double const* const p1(nodes[1]->getData());
        double const* const p2(nodes[2]->getData());
        const double v1[3] = {p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]};
        const double v2[3] = {p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2]};
        CrossProduction(v1, v2, normal_vector);
        NormalizeVector(normal_vector, 3);
    }
    if (this->GetElementType() ==
        MshElemType::QUAD)  // BG 04/2011: calculation of the normal vector of a
                            // quad element
    {
        double const* const p0(nodes[0]->getData());
        double const* const p1(nodes[1]->getData());
        //		double const* const p2 (nodes[2]->getData()); // TF unused
        // variable
        double const* const p3(nodes[3]->getData());
        const double v1[3] = {p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]};
        const double v2[3] = {p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2]};
        // v1[0] = nodes[1]->X() - nodes[0]->X();
        // v1[1] = nodes[1]->Y() - nodes[0]->Y();
        // v1[2] = nodes[1]->Z() - nodes[0]->Z();
        // v2[0] = nodes[3]->X() - nodes[0]->X();
        // v2[1] = nodes[3]->Y() - nodes[0]->Y();
        // v2[2] = nodes[3]->Z() - nodes[0]->Z();
        CrossProduction(v1, v2, normal_vector);
        NormalizeVector(normal_vector, 3);
    }
}

// KR 2010/11/16
void CElem::setElementProperties(MshElemType::type t)
{
    switch (t)
    {
        case MshElemType::LINE:
            nnodes = 2;
            nnodesHQ = 3;
            ele_dim = 1;
            geo_type = MshElemType::LINE;
            nfaces = 2;
            nedges = 1;
            break;
        case MshElemType::QUAD8:
            nnodes = 4;
            nnodesHQ = 8;
            ele_dim = 2;
            geo_type = MshElemType::QUAD8;
            nfaces = 4;
            nedges = 4;
            break;
        case MshElemType::QUAD:
            nnodes = 4;
            nnodesHQ = 9;
            ele_dim = 2;
            geo_type = MshElemType::QUAD;
            nfaces = 4;
            nedges = 4;
            break;
        case MshElemType::HEXAHEDRON:
            nnodes = 8;
            nnodesHQ = 20;
            ele_dim = 3;
            nfaces = 6;
            nedges = 12;
            geo_type = MshElemType::HEXAHEDRON;
            break;
        case MshElemType::TRIANGLE:
            nnodes = 3;
            nnodesHQ = 6;
            ele_dim = 2;
            geo_type = MshElemType::TRIANGLE;
            nfaces = 3;
            nedges = 3;
            break;
        case MshElemType::TETRAHEDRON:
            nnodes = 4;
            nnodesHQ = 10;
            ele_dim = 3;
            geo_type = MshElemType::TETRAHEDRON;
            nfaces = 4;
            nedges = 6;
            break;
        case MshElemType::PRISM:
            nnodes = 6;
            nnodesHQ = 15;
            ele_dim = 3;
            geo_type = MshElemType::PRISM;
            nfaces = 5;
            nedges = 9;
            break;
        case MshElemType::PYRAMID:
            nnodes = 5;
            nnodesHQ = 13;
            ele_dim = 3;
            geo_type = MshElemType::PYRAMID;
            nfaces = 5;
            nedges = 8;
            break;
        default:
            std::cerr << "CElem::setElementProperties MshElemType not handled"
                      << "\n";
    }
    this->nodes_index.resize(quadratic ? nnodesHQ : nnodes);
}

// NW
double* CElem::ComputeGravityCenter()
{
    const size_t nnodes0 = this->nnodes;
    for (size_t i = 0; i < nnodes0; i++)  // Nodes
    {
        double const* const pnt(nodes[i]->getData());
        gravity_center[0] += pnt[0];
        gravity_center[1] += pnt[1];
        gravity_center[2] += pnt[2];
    }
    gravity_center[0] /= (double)nnodes0;
    gravity_center[1] /= (double)nnodes0;
    gravity_center[2] /= (double)nnodes0;
    return gravity_center;
}

/**************************************************************************
MSHLib-Method:
11/2014 JOD Implementation
**************************************************************************/
void CElem::DirectNormalVector()
{
    for (int i = 0; i < 3; i++)
        if (normal_vector[i] < 0)
            InvertNormalVector();
}

/**************************************************************************
MSHLib-Method:
11/2014 JOD Implementation
**************************************************************************/
void CElem::InvertNormalVector()
{
    normal_vector[0] = -normal_vector[0];
    normal_vector[1] = -normal_vector[1];
    normal_vector[2] = -normal_vector[2];
}

}  // namespace MeshLib
